clear all

% Simulate benchmark levels

global Dividend

Dividend=0;
%
% Given Par2 optimal:
n = 50; tau = 0.1;
r = 0.036;
Par1 = zeros(25,1);
load ParEst_eh0_time0_Mex ParEst_eh0_time0_Mex
load ParEst_eh0_time1_Mex ParEst_eh0_time1_Mex
Par1(1:24) = ParEst_eh0_time0_Mex;
Par1(21) = 10^Par1(21);
Par1(22:24) = Par1(22:24)*1000;
Par1(25) = 0;
gamma =  Par1(25); a = Par1(22); alphaf = Par1(17); betaf = Par1(18); alphai = Par1(19); betai = Par1(20); theta = Par1(21);
b1 = Par1(23); b2 = Par1(24);
load W_time1_age0_eh0_comp11.out
w_f_g=W_time1_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn,Rw_fi_ni, Rw_ii_ni, Rw_ii_in, Rw_if_in,...
     Rw_fi_fn, Rw_ff_fn, Rw_if_nf, Rw_ff_nf, Rw_in_nn, Rw_fn_nn, Rw_ni_nn, Rw_nf_nn]=solve_value_function_A3_C(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
[HHstates,wage_moments,transitions,gff, gfi, gif, gii, gfn, gnf, gin, gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn);
wage_moments_gamma_0_GE = wage_moments;
wage_moments_gamma_0_GE_red=wage_moments(1:8);%meanwi_1; meanwf_1; meanwi_2 ; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2

HH_states_gamma_0_GE = HHstates;
save wage_moments_gamma_0_GE_red wage_moments_gamma_0_GE_red
save HH_states_gamma_0_GE HH_states_gamma_0_GE

mff=HHstates(1);
mfi=HHstates(2);
mfn=HHstates(3);
mif=HHstates(4);
mnf=HHstates(5);
mii=HHstates(6);
min_=HHstates(7);
mni=HHstates(8);
mnn=HHstates(9);
Un_Head = mnf + mni + mnn; Un_Spouse = mnn + min_ + mfn; Inf_rate = (mnn + min_ + mni + mii);
Agg_stocks_gamma_0_GE = [Inf_rate; Un_Head; Un_Spouse];
gf_1 = zeros (50,1); gf_2 = zeros (50,1); gi_1 =  zeros (50,1); gi_2 = zeros (50,1);
 for i = 1:50
 gf_1(i) =  (mff*(sum(gff(i,:))) + mfi*(sum(gfi(i,:))) + mfn*gfn(i))/(mff+mfi+mfn);   
 end

  for i = 1:50
 gi_1(i) =  (mif*(sum(gif(i,:))) + mii*(sum(gii(i,:))) + min_*gin(i))/(mif+mii+min_);   
  end
 
   for i = 1:50
 gf_2(i) =  (mff*(sum(gff(:,i))) + mif*(sum(gif(:,i))) + mnf*gnf(i))/(mff+mif+mnf);   
   end
 
    for i = 1:50
 gi_2(i) =  (mfi*(sum(gfi(:,i))) + mii*(sum(gii(:,i))) + mni*gni(i))/(mfi+mii + mni);   
 end
save Agg_stocks_gamma_0_GE Agg_stocks_gamma_0_GE

% Transitions 
Trans_head = [transitions(3); transitions(4); transitions(1); transitions(5); transitions(2); transitions(6); transitions(13); transitions(14)];
Trans_spouse = [transitions(9); transitions(10); transitions(7); transitions(11); transitions(8); transitions(12); transitions(15); transitions(16)];
Transitions_gamma_0_GE = [Trans_head;NaN;NaN;Trans_spouse];
save Transitions_gamma_0_GE Transitions_gamma_0_GE



